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ABSTRACT 

In this article, we present a novel scheme for segmenting the im¬ 
age boundary (with the background) in optoacoustic small animal 
in vivo imaging systems. The method utilizes a multiscale edge de¬ 
tection algorithm to generate a binary edge map. A scale depen¬ 
dent morphological operation is employed to clean spurious edges. 
Thereafter, an ellipse is fitted to the edge map through constrained 
parametric transformations and iterative goodness of fit calculations. 
The method delimits the tissue edges through the curve fitting model, 
which has shown high levels of accuracy. Thus, this method enables 
segmentation of optoacoutic images with minimal human interven¬ 
tion, by eliminating need of scale selection for multiscale processing 
and seed point determination for contour mapping. 

Index Terms — Image segmentation. Photoacoustic effects. Im¬ 
age edge detection, Curve fitting 

1. INTRODUCTION 

Optoacoustic Tomography (OAT) has emerged in the last decade 
as a new hybrid bioimaging modality combining the advantages of 
acoustical detection and optical absorption contrast. OAT provides 
structural as well as functional information of biological tissues in 
two (2D) and three (3D) dimensions, with resolution and frame rates 
representative of ultrasound (T). Research efforts in OAT have been 
directed towards the development of new hardware components and 
inversion methodologies allowing increasing imaging speed and res¬ 
olution, as well as on investigating potential biomedical applications. 
Recently developed OAT cross-sectional systems enable whole-body 
small animal in-vivo imaging. The unique capabilities offered open 
up the unexplored domain of post-reconstruction image analysis for 
this modality. Segmentation is one the most challenging tasks in 
image processing, and the relatively low intrinsic contrast of back¬ 
ground structures (compared to digital photography) and limited 
view problems of OAT increase the complexity involved. Classical 
experiments on non-human primates demonstrate that complex ob¬ 
jects are perceived by vision system in a multiscale manner, and are 
tuned from coarse scale to finer scales m , we use the information as 
a basic assumption in our study. For the edge detection at individual 
scales, the Sobel operator, which approximates the gradient of the 
image intensity function, or the Canny edge detector El which uses 
a feature synthesis step from fine to coarse scales, are commonly 
employed. Multiscale edge detection, a rich area of interest in itself, 
can enhance the performance of edge detection in a fixed (original) 
scale. For example, the Perona-Malik flow (anisotropic diffusion) |4| 
gives a successful formulation of scale-space adaptive smoothing 
function, leading to better preservation of edges while achieving a 


better noise performance. These edge detection methods have been 
demonstrated to be useful for optoacoustic images and employed for 
calibration of reconstruction parameters by Mandal et.al.(5)- More 
sophisticated techniques were proposed by Tabb and Ahuja and 
followed by Ma and Manjunath (Tjwho developed the concept of 
designing a vector field for edge detection. In El the vector field 
is created by analyzing the neighborhood of a pixel, where the size 
and spatial scale of the neighborhood were determined heuristically 
by an homogeneity parameter, and is adaptively determined on a 
pixel to pixel basis. The edge flow method suggested by |f7| utilizes 
the color and textural information of images to track changes in 
directions, creating a vector flow. This method detects boundaries 
when there are two opposite directions of flow at a given location 
in a stable state. However, it depends strongly on color information 
and requires a user defined scale to be input as a control parameter. 
In a follow up work by El > a methods was demonstrated which 
eliminates the need of scale selection and potentially works with 
gray-scale images. We have developed a new segmentation method 
for OAT by integrating the multiscale edgeflow, scale space diffu¬ 
sion and morphological image processing. The suggested multiscale 
edge detection method is illustrated in section 2, and the processing 
results of OAT images are displayed in section 3. In section 2.4, 
we discuss about delineating the tissue boundaries in OAT, which 
is another challenging issue we aim to tackle. Due to illumina¬ 
tion conditions and limited view considerations in the detection, 
the reconstructed images can be affected by open edges and fuzzy 
boundaries. Traditional edge linking algorithms may then fail to 
provide satisfactory segmentation. Thereby, we suggest a novel 
curve fitting model that delimits the tissue-background boundaries 
in OAT images. We further characterize the goodness of fit (GoF) 
by using Dice coefficient metric (DM), an standard evaluation ma¬ 
trix. The GoF is employed to iteratively improve the segmentation 
performance. 

2. METHODS AND ALGORITHMS 

The proposed methodology aims at providing a framework of OAT 
image segmentation and reducing the parameters that need to be de¬ 
fined to achieve a segmented boundary between imaged biological 
tissue and acoustical coupling medium. Details on the imaging setup 
and the implementation steps of the algorithm are given in the fol¬ 
lowing subsections. 

2.1. Imaging Setup and Protocol 

A commercial whole-body small animal multispectral optoacoustic 
tomography (MSOT) scanner (Model MSOT256-TF, iTheraMedical 


GmbH, Munich, Germany) was used for data acquisition. The scan¬ 
ner consists of a 256-element array of cylindrically focused piezo¬ 
composite transducers with 5 MHz central frequency. The trans¬ 
ducer array covers an angle of approximately 270 degrees and has 
a radius of curvature of 40 mm. Light excitation is provided with 
the output laser beam from a wavelength-tunable optical parametric 
oscillator (OPO)-based laser, which is shaped to attain ring-type uni¬ 
form illumination on the surface of the animal (mouse) by means of 
a fiber bundle with an output arranged on annular illumination slots. 
The system allows simultaneous acquisition of the signals generated 
with each laser pulse, which are digitized at 40 MS/s. Thereby, the 
scanner is capable of rendering 10 cross-sectional images per sec¬ 
ond. For acquisition of the entire mouse cross-sections the laser 
wavelength was set between 690 and 900 nm and the signals were 
averaged 10 times in order to improve SNR performance. The op- 
toacoustic images representing the distribution of optical absorption 
were reconstructed using a non-negative constrained model-based 
inversion algorithm (8). Prior to reconstruction, the signals were ini¬ 
tially band-pass filtered with cut-off frequencies between 0.1 and 7 
MHz for removing low-frequency offsets and high-frequency noise. 
The reconstructed images have a 200 x 200 effective pixel resolu¬ 
tion (for a reconstruction area of 20 x 20 mm^) and are represented 
in gray-scale (normalized 0-255 level for image processing appli¬ 
cations). The animal handling for imaging applications were per¬ 
formed in full conformity with institutional guidelines and with ap¬ 
proval from the Government of Upper Bavaria (Germany). 

2.2. Edge-flow vector fleld and edge detection 

The general trend in multiscale edge detection is to define a scale 
a-priori and then estimate the scale locally, however, Sumengen and 
Manjunath O suggested a geometrically inspired method that esti¬ 
mates the edges that exists both in coarse and fine scales, and localize 
them in the fine scale. In this article, we propose a modified version 
of the edgeflow methods GHH, by integrating anisotropic diffusion 
and scale-space dependent morphological processing for convenient 
application on OAT images. The edge flow algorithm defines a vec¬ 
tor field, such that the vector flow is always directed towards the 
boundary on both its sides. The classical edge flow model utilizes a 
vector propagation stage. In the current method, the relative direc¬ 
tional differences are considered for computing gradient vector. The 
gradient vector strengthens the edge locations and tracks the direc¬ 
tion of the flow along x and y directions. The search function looks 
for sharp changes from positive to negative signs of flow directions 
and whenever it encounters such changes, the pixel is labeled as an 
edge point. The magnitude of the change is the deciding factor be¬ 
hind the edge strength, which is reflected as edge intensity in the 
final edge map. The vector field is generated explicitly from fine to 
coarse scales, whereas the multiscale vector conduction is implicitly 
from coarse to finer scales. Thus, the algorithm is suitable for local¬ 
izing the edges in the finer scales, which is achieved by preserving 
only the edges and neighborhoods that exists in several scales (de¬ 
pending on the threshold employed), and suppressing features that 
disappear rapidly with increment of scales. The workflows for the 
vector field generation and morphological processing are illustrated 
in the Algorithm 1. The pseudocode further outlines the process of 
geometric curve fitting. In the current implementation, the Gaussian 
offset (cr) is also lowered for the gray-scale OAT images. For the 
acquired in vivo OAT images, optimum performance was achieved 
with (T = 3. We analyze the images between scales s=l, and s=3, 
where s=l is the starting (finest) scale. The interval is sampled at 
sub-pixel resolution [As = 0.5] for tracing the dislocation of edges 


in subsequent edges, as used by Berghold (9). 


Algorithm 1 Multiscale Edge Detection and Curve Fit 
I(x,y) <r- image 

si, Sn <r- smallest and largest scale of the image respectively 
As = 0.5 be the sampling interval for each scale 
s = si > initialize book keeping 

U, Unew ^ VI{x,y) > Gradient Vector 

Une.w Vc.'VI{x,y) + c{x,y,t)AI t> Anisotropic Diffusion 
while Si < Sn do 
s •<— s -F As 
M <— maa;(|| tJ ||) 

Unew •<— V/(a;, y) at scale s 
Unew •<- Vc.VI{x,y) + c{x,y,t)^I 
for <each pixel in image> do 

if II U{x,y) II < M/C then o C is thresholding constant 
U{x,y) = Unew{x,y) 

elseii abs{arctan{U{x,y),Unew{x,y))) < vr/dthen 
U {X, y) = U {X, y) + Unew (X, t/) 

else 

U{x,y) = U{x,y) 

end {if, for} 

The final U is the edge flow vector we are interested in 
U <— Binarize(U) 

Enhance U by performing morphological operations 
for gradient magnitude image in scale Sn do 
if n > 2, use strel— > Opx, disc 

else t> structural element)strel) 

strel— > Ipx, disc then 

Apply Erosion operation; Close to recover edges 

end for 

Obtain centroids from binary edge-map 

Obtain minBoundCircle and minBoundEllipse on centroids 

Calculate Dice Coefficient > Iterate over scales 

end while 

The final boundary is given by minBoundCircle and or minBound¬ 
Ellipse after GoE maximization. 


The algorithm searches the edges in finer scales and strengthen 
them with the edges recovered from higher scale. The homogeneous 
regions have vectors of zero length, so the detected edge segments 
grow in thickness (and often strength) as we move from lower to 
higher scales. Some edges do not exist in lower scale, but can still 
be significant. To decide on the same and reduce noise we put a 
boundary condition - and when the maximal edge strength (M) is 
greater than a heuristically predefined constant (C) -the edges are re¬ 
tained, or else they are discarded.The primary objective of applying a 
edge detection is to delineate the boundary of the imaged object and 
differentiate it from the background. But often in optoacoustics, the 
signal originate from the impurities or inhomogenities withing cou¬ 
pling medium. Further, noisy background is present in reconstructed 
images (lower boundary in Fig l.a) due to limited view, and short¬ 
comings of inversion methodologies. This noises are often strong 
enough to be detected by edge detection algorithm as true edges. 
Thus, we use an anisotropic diffusion process to further clean up 
the image, the diffusion process smoothens the image without sup¬ 
pressing the edges. Thereafter a non-linear morphological process¬ 
ing is done on the binary (diffused)edge mask. We take an sub-pixel 
sampling approach (0.5 px), rendering the operation is redundant 
beyond the second scale level. The morphological mask is differen- 





tially chosen at different scales. Initially, the image is eroded with 
a disc structuring element to remove noisy patches, but it also thins 
the edges. To recover the edges a closing operation is executed, with 
smaller structural element for erosion and a bigger element (2px) for 
dilation. 


2.3. Parametric shape model fitting and Goodness of fit 

In sectiot i2^ we discussed about detection of edges, several authors 
have utilized the edge-flow vector for segmentation. In OAT im¬ 
ages we see formation of smaller edge clusters and open contours. 
Thus, getting a ideal segmentation using edge linker seem to perform 
poorly. However, given the fact that our current problem which re¬ 
quires segmenting the image into only two classes - image and back¬ 
ground, we follow a simple curve fitting approach. The proposed 
method first generates the centriods for edge clusters and then try 
to fit on a geometric pattern (deformable ellipse) iteratively through 
a set of parametric operations. A typical scenario in curve fitting 
is when the data is a best fitted, but some data points lies outside 
the curve, as we take an interpolated spline fitting criterion. In our 
approach, we modeled an inclusion criterion which enclosed all cen¬ 
troids and create forms a closed curve. Theoretically, it draws a 
convex hull with the centroids on a perimeter and approximate it 
to the nearest curve. In some datasets, we see presence of the out¬ 
lier centroids which significantly biases the curve (due to presence 
of the inclusion criterion) leading to erroneous fitting. To avoid such 
complications, the values of the centroid positions are filtered for 
squeezing out outliers through a median filter, and then analyzed for 
generating the shape models which match the object boundary.The 
goodness of fit (GoF) is calculated using the average DM, which is 
a measure of contour overlap utilizing the area under the fitted curve 
(A), manually segmented region (M), and their intersection. DM al¬ 
ways vary between 0 and 1, with DM > 7 being considered a good 
segmentation (To). The DM can be expressed as; 


DM 


2 I (A n M) I 
IA|+|M| 


( 1 ) 


The GoF measure is used to iteratively improve the segmentation 
performance - the images are investigated at deeper scale spaces and 
the morphological diffusion step is decided based on the GoF value. 
We break away from the loop once a stable GoF is obtained and 
going to lower scale space negatively impact the calculated GoF val¬ 
ues.The values for Rand Index (RI), which gives a measure of sim¬ 
ilarity and statistically used to quantify accuracy are also calculated 
and plotted in Figure 3. 


3. RESULTS AND DISCUSSION 

For testing and standardization of the introduced image processing 
framework, a database of 30 datasets for 3 different anatomical re¬ 
gions in mice (10 for each region) was created. Each individual 
dataset represents signals (averaged 10 times) acquired at up to 6 
different positions and at 8 different wavelengths (between 690-900 
nm). Since manual segmentation of each dataset was needed for 
computing the performance of the algorithm vis-a-vis ground truth, 
the number of datasets was scaled down. Thereby, 4 datasets were 
considered from each anatomical regions, namely brain, liver and 
kidney/spleen, for computing the goodness of fit parameters. The 
proposed algorithm demonstrated good edge recovery performance. 
As shown in Fig 1, the algorithm weighted the edges that appear 
in multiple scales and reduced the spurious edges. The assumption 
made by (7l| is that if the vector direction change matches in multiple 


scales, then we can infer that it corresponds to real edges. Thus, the 
vector directions were checked in both finer and coarser scales, in a 
way that when there is a match, the designated edges were strength¬ 
ened. The construction of the algorithm allows to detect the initial 
edge points from the finer scales, and reinforce then as we move to 
the larger scales. The finer scales are more immune to noise and the 
use of a non-negative constraint during the image reconstruction pro¬ 
cess prevents unnatural movement of the vector field (potentially due 
to absence of undesired negative values). Thus the edges detected in 
the finer scales are very significant to recover object boundaries, and 
are helpful in segmentation of OAT images. In Figure 2, we show 
the performance of multiscale segmentation (2c) along with the edge 
map recovered using Sobel operator (2b). The improvement in the 
edge detection performance by considering multiple scale (as over 
single scale in Sobel) is evident in the combined multiscale edge 
map ((7 = 1 — 3). Further, a closer observation reveals that the 
morphological processing have successfully eliminated the spurious 
edges formed beyond the tissue boundary (Fig. 2d). Finally, in Fig¬ 
ure 2(d-e) we show the calculated centroid clusters obtained from the 
morphologically processed binary edge map, and the ellipse fitting 
model applied to this centroids respectively. 

Thereafter, we computed the goodness of fit using quantitative 
measures, viz. DM and RI, with ground truth (manual segmen¬ 
tation) as reference. In Table 1 we show the performances of the 
curve fitting using the DM 0. Theoretically, DM values above 
0.7 are considered to represent a good segmentation result, and us¬ 
ing the proposed methodology DM values between 0.90 and 0.96 
were achieved. Morphological processing and increasing scale- 
space depth was employed by tracking the corresponding changes 
in GoF value to improve noise performance when outlier centroids 
are present. Improvements in DM values was observed in most 
regions by adding anisotropic diffusion and morphological filtering 
with iterative GoF optimization, except for the kidney/spleen zone 
where the DM decreased after filtering, although DM values above 
0.95 both with and without secondary filtering were achieved. From 
figure 3, it is be clearly inferred from both DM and RI that the 
improvements in the curve fit performance is observed in scale 2 
followed by a marginal change in scale 3. So, for most practical pur¬ 
poses and reducing computational overhead, as scale depth of 2 can 
be chosen for OAT imaging experiments and image segmentation. 



Fig. I. Reconstructed optoacoustic image using non-negative con¬ 
strained model based inversion (a), and (b-d) edges detected at mul¬ 
tiple scales s = 1, 2, 3 respectively. 





Fig. 2. Reconstructed Optoacoustic image with non-negative con- 
strain(a), edge-maps generated by Sobel operator (b), Multiscale 
(Edgeflow) map (c)and (proposed) Morphological processed multi¬ 
scale edge map (d) are shown. In (e) computed centroids from edge 
clusters are shown in blue, and (f) illustrates the curve fitting method 
applied to the centroids (fitted curve marked in red). [Scale = 3mm] 


Table 1. Efficacy of Curve Fitting 


Regions 

Dice Coefficient Metrics 
Non-optim GoF-Optim 

Rand Index 

Brain 

0.9387 

0.9440 

0.9771 

Liver 

0.9093 

0.9447 

0.9396 

Kidney/Spleen 

0.9606 

0.9571 

0.9672 


4. CONCLUSIONS 

In this article, a method of delineating the boundary of optoacoustic 
small animals images using a multiscale edge detection algorithm in 
combination with geometrical curve fitting has been presented. It is 
noteworthy that the method employed is self-deterministic and re¬ 
quires minimal human intervention. The optoacoustic signal/image 
datasets can be very large given the 5D nature of this modality (m, 
automating the image formation and analysis workflows is a very 
challenging and important problem. Thus, the algorithms and the 
workflows demonstrated herein are expected to be helpful in au¬ 
tomating optoacoustic image segmentation, with important signifi¬ 
cance towards enabling quantitative imaging applications. 
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Eig. 3. GoF trace and the standard deviation for the 
torso(kidney/spleen region) scans of mouse. (Total dataset size = 
12 (z-slices in whole body small animal tomographic scanner), ac¬ 
quired from 2 different mice) in-vivo 
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